---
title: "Black Sea sprat XSA stock assessment (`r as.numeric(data.stock@range['minyear'])`-`r as.numeric(data.stock@range['maxyear'])`)"
author: "Author: Piatinskii M., Azov-black sea branch of VNIRO"
date: 'Report build date: `r Sys.time()`'
output:
  html_document:
    toc: true
    toc_depth: 3
    toc_float: true
    theme: default
---

```{r echo=FALSE}
year.terminal <- as.numeric(data.stock@range["maxyear"])
year.start <- as.numeric(data.stock@range["minyear"])
```

## 1. Input data 
Here you can review input data summary information. This input data used to perform XSA model and make short-term forecasting.

### 1.1 Stock (FLStock)

```{r echo=TRUE}
summary(data.stock)
```

### 1.2 Index (FLIndices)

```{r}
summary(data.indices)
```

### 1.3 Parametrization

```{r echo=FALSE}
for (var in ls())
  if (startsWith(var, "config.") & var != "config.tuning.xsa") print(sprintf("%s=%s", var, environment()[[var]]))
```

## 2. Preliminary diagnostic
In this section you should review all available diagnostic information before start analysis. If diagnostics show big residuals, different slope of log index regression, retrospective miss convergence - you should **stop** analysis with "bad quality" input data assumption. The procedure should include next steps:

1. Review shrinkage impact factor: how shrinkage Fse change results?
2. Review retrospective stability at diffenrent shrinkage Fse levels
3. Review Mohn's-Rho retrospective test results
4. Review Survey age-vs-age regression lines: all lines should have the same slope (b coef) to be representative by time vector

### 2.1 Shrinkage impact factor
Shrinkage standard error (Fse, or shrinkage s.e.) is a major model covergence factor for model fitting to observed data. By default, Fse = 0.5. 

```{r fig.cap="Fig.2.1.1. Full XSA estimates at different Fse levels", echo=FALSE}
knitr::include_graphics("output/1_FseChoose/summary_all.png")
```

```{r fig.cap="Fig.2.1.2. Last 5 years XSA estimates at different Fse levels", echo=FALSE}
knitr::include_graphics("output/1_FseChoose/summary_short.png")
```

```{r echo=FALSE}
d <- data.frame(year = (year.start):(year.terminal), 
                fse0.5 = as.vector(ssb(result.fse.summary[[1]])),
                fse1.0 = as.vector(ssb(result.fse.summary[[2]])), 
                fse1.5 = as.vector(ssb(result.fse.summary[[3]])), 
                fse2.0 = as.vector(ssb(result.fse.summary[[4]])), 
                fse2.5 = as.vector(ssb(result.fse.summary[[5]]))
                )

d2 <- data.frame(year = (year.start):(year.terminal), 
                fse0.5 = as.vector(fbar(result.fse.summary[[1]])),
                fse1.0 = as.vector(fbar(result.fse.summary[[2]])), 
                fse1.5 = as.vector(fbar(result.fse.summary[[3]])), 
                fse2.0 = as.vector(fbar(result.fse.summary[[4]])), 
                fse2.5 = as.vector(fbar(result.fse.summary[[5]]))
                )

knitr::kable(d, caption = "Table 2.1.1. SSB estimates at different Fse levels")
knitr::kable(d2, caption = "Table 2.1.2. Fbar estimates at different Fse levels")
```

### 2.2 Retrospective stability
There is classic restrospective analysis with step-by-step reducing input time series length by 1 year. The main target of retrospective analysis is observe model convergence in luck-of-data terms. 

```{r fig.cap="Fig.2.2.1. Retrospective at Fse = 0.5", echo=FALSE}
knitr::include_graphics("output/2_RetroFse/0.5.png")
```

```{r fig.cap="Fig.2.2.2. Retrospective at Fse = 1.0", echo=FALSE}
knitr::include_graphics("output/2_RetroFse/1.png")
```

```{r fig.cap="Fig.2.2.3. Retrospective at Fse = 1.5", echo=FALSE}
knitr::include_graphics("output/2_RetroFse/1.5.png")
```


```{r fig.cap="Fig.2.2.4. Retrospective at Fse = 2.0", echo=FALSE}
knitr::include_graphics("output/2_RetroFse/2.png")
```

```{r fig.cap="Fig.2.2.5. Retrospective at Fse = 2.5", echo=FALSE}
knitr::include_graphics("output/2_RetroFse/2.5.png")
```

### 2.3. Mohn's-Rho test
The basic ICES (2018) procedure to determine model stability over years - retrospective Mohn rho test. Mohn's rho test calculate relative bias for latest 3 years (in default ICES procedure - 5 years) retrospective variations on scale -1 ... +1. Low negative values of *rho* leads to underestimate factor, high positive values - to overestimation (remark: values lowest -0.4 or higher +0.4 shows high variations and low model stability). So procedure is:

$$
\begin{aligned}
relbias = (retro - base) / base \\
rho = mean(relbias)
\end{aligned}
$$

For SSB and F math approach will be

$$
\begin{aligned}
\rho_{SSB} = \frac{1}{n} \sum_{i=-5}^{n=0} \frac{(SSB_{i} - \overline{SSB})}{\overline{SSB}} \\
\rho_{Fbar} = \frac{1}{n} \sum_{i=-5}^{n=0} \frac{(Fbar_{i} - \overline{Fbar})}{\overline{Fbar}} \\
\end{aligned}
$$

where i - year steps from terminal year (terminal-3, terminal-2, ... terminal).

```{r echo=FALSE}
result.retro.diagnostic.ssb %>%
  rownames_to_column(., var="year") %>%
  mutate_if(is.numeric, ~round(.,1)) %>%
  knitr::kable(., caption = "Table 2.3.1. SSB retrospective values", row.names = TRUE)
```

```{r echo=FALSE}
result.retro.diagnostic.f %>%
  rownames_to_column(., var="year") %>%
  mutate_if(is.numeric, ~round(.,3)) %>%
  knitr::kable(., caption = "Table 2.3.2. Fbar retrospective values", row.names = TRUE)
```

**Final** Rho estimates: 

$$\rho_{SSB} = `r round(result.retro.mohnrho.ssb,3)`$$
$$\rho_{Fbar} = `r round(result.retro.mohnrho.fbar,3)`$$


### 2.4 Survey index regressions
Surveys index regression diagnostics. The same regressions slope show a good quality of survey tuning. Survey vs survey plot (fig. 2.4.3) indicate fish numbers distribution between surveys. High correlation r2 value leads to good survey compatibility.

```{r fig.cap="Fig.2.4.1. Survey 1 log index in time vector", echo=FALSE}
knitr::include_graphics("output/4_SurveysRegressions/survey_1.png")
```

```{r fig.cap="Fig.2.4.2. Survey 2 log index in time vector", echo=FALSE}
knitr::include_graphics("output/4_SurveysRegressions/survey_2.png")
```

```{r fig.cap="Fig.2.4.3. Survey vs survey age-to-age comparision", echo=FALSE}
knitr::include_graphics("output/4_SurveysRegressions/summary_all.png")
```


## 3. Model final tuning (FLXSA.control)
There is a final FLXSA.control object with tuning used in XSA fitting.

```{r echo=FALSE}
print(config.tuning.xsa)
```

## 4. XSA results
Here the final short result of XSA model fitting. More detailed information can be found in appendix.

### 4.1. Basic results

```{r echo=FALSE}
df <- data.frame(year = year.start:year.terminal, 
                 ssb.tons = as.vector(ssb(result.xsa.stock)),
                 rec = as.vector(rec(result.xsa.stock)),
                 fbar = round(as.vector(fbar(result.xsa.stock)), 3)
                 )

knitr::kable(df, caption = "Table 4.1. XSA estimate results")

```

Survivour numbers at age distribution is a major estimate in XSA analysis. Based on stock numbers and fishing mortality information XSA calculate SSB, Fbar and Recruitment numbers. Fig. 4.1.1 show the stock numbers at age estimation. Fig 4.1.2 show the fishing mortality at age estimates. Fig. 4.1.3 show calculated summary results based on stock numbers: spawning stock biomass(SSB), recruitment numbers (rec) and Fbar estimates.

```{r fig.cap="Fig.4.1.1. Stock numbers at age by time vector", echo=FALSE}
knitr::include_graphics("output/5_XsaEstimates/stock_numbers_at_age.png")
```

```{r fig.cap="Fig.4.1.2. Fishing mortality at age by time vector", echo=FALSE}
knitr::include_graphics("output/5_XsaEstimates/f_at_age.png")
```

```{r fig.cap="Fig.4.1.3. Summary XSA results", echo=FALSE}
knitr::include_graphics("output/5_XsaEstimates/summary_all.png")
```

### 4.2. Uncertainty variations
In this section XSA results presented under uncertainty investigation. The uncertainty has been identified based on FAO approach as st.dev of fishing mortality destribution. The uncertainty is: 
$$
\begin{aligned}
\sigma = \sqrt{\frac{1}{N-1} \sum_{i=1}^N (F_i - \overline{F})^2}   \; [1]  \\
\end{aligned}
$$
Uncertainty justification (sigma) used to determine variations in ssb, rec and F as linear multiplier:
$$
\begin{aligned}
var[SSB]_{i} = SSB_{i} *  \sigma  \; [2]  \\
var[REC]_{i} = REC_{i} *  \sigma  \; [3]  \\
var[F]_{i} = F_{i} *  \sigma  \; [4]  \\
\end{aligned}
$$
Then variation used to create norman distribution of estimates based on mean and sd. The final results shown at fig. 4.2.1 

```{r fig.cap="Fig. 4.2.1. Uncerteinty factor for Rec, SSB, F", echo = FALSE}
knitr::include_graphics("output/5_XsaEstimates/summary_uncertainty.png")
```

In plot technique the classic quantile method used for 0.05(lower), 0.5 (mean) and 0.95 (upper) levels.

## 5. XSA diagnostics
In process of XSA model fit iterations survey information used to tune stock data. Observed vs fitted values (with survey) should be reviewed as fitting quality factor. High residual values leads to bias in model. In fact values in range -2.25 < log(residual) < 2.25 are not significant. Values out of this range idicate low survey quality (by single year or age group). Full negative (year or age group) residual indicate high posibility of bias. 

```{r fig.cap="Fig.5.1. Log residual diagnostics for survey1 vs model fit", echo=FALSE}
knitr::include_graphics("output/6_Residuals/summary_survey1.png")
```

```{r fig.cap="Fig.5.2. Log residual diagnostics for survey2 vs model fit", echo=FALSE}
knitr::include_graphics("output/6_Residuals/summary_survey2.png")
```

```{r fig.cap="Fig.5.3. Log catchability residuals at age diagnostics", echo=FALSE}
knitr::include_graphics("output/6_Residuals/catch_at_age_summary.png")
```

Log-residual values raw output:
```{r echo=FALSE}
print(result.xsa.fit@index.res)
```

Min/max residuals:
```{r echo=FALSE}
for (sur in names(result.xsa.fit@index.res)) {
  print(paste("Survey: ", sur))
  cat(paste("Min log.res =", min(result.xsa.fit@index.res[[sur]]), "\n"))
  cat(paste("Max log.res =", max(result.xsa.fit@index.res[[sur]]), "\n"))
}
```